Efficient simulation of strong system-environment interactions 
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Multi-component quantum systems in strong interaction with their environment are receiving increasing at- 
tention due to their importance in a variety of contexts, ranging from solid state quantum information processing 
to the quantum dynamics of bio-molecular aggregates. Unfortunately, these systems are difficult to simulate as 
the system-bath interactions cannot be treated perturbatively and standard approaches are invalid or inefficient. 
Here we combine the time dependent density matrix renormalization group methods with techniques from the 
theory of orthogonal polynomials to provide an efficient method for simulating open quantum systems, includ- 
ing spin-boson models and their generalisations to multi-component systems. 
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Introduction- Dissipation and decoherence caused by 
noisy, fluctuating environments are fundamental and ubiq- 
uitous phenomena which appear in almost all experiments 
on quantum systems. While initially thought to lead to the 
rapid erasure of quantum signatures, it was later realized that 
the presence of environmental noise can in fact be instru- 
mental for maintaining quantum correlations in the steady 
state U]. Decoherence is also thought to be essential for ex- 
plaining the efficient excitation energy transport (EET) in the 
pigment-protein complexes (PPCs) of photosynthetic organ- 
isms 12l|3|], biological systems in which long-lived quantum 
coherence has been recently observed J3]. As a result, novel 
emphasis is being placed in better understanding the inter- 
play between the coherent quantum dynamics and incoher- 
ent processes arising from the interaction with the environ- 
ment yfl. Ultimately, this understanding may facilitate the 
development of noise-assisted devices and, potentially, light- 
harvesting systems with increased efficiency. 

To date, the dynamical behaviour of interacting open quan- 
tum systems and, in particular, noise-assisted transport, is 
frequently investigated in terms of simple dynamical models 
in which environmental dephasing and relaxation are treated 
with Lindblad or Bloch-Redfield master equations. These 
standard methods are both based on the assumptions of weak 
system-bath coupling and the Markov approximation. How- 
ever, these approximations are not valid for many realistic sys- 
tems of current interest, and assuming that the correlation time 
of the environments in these examples is much faster than the 
system dynamics is frequently not justified. For instance, in 
typical PPCs the dynamical timescales of the bath can be com- 
parable or even slower than the EET dynamics 10,0]. More- 
over, in the limit of slow bath dynamics, perturbative treat- 
ments of the system-environment coupling cannot be used 
even if the system-bath coupling is intrinsically weak J^,§]. 
Recently, important steps have been taken towards develop- 
ing and analysing fully non-perturbative and non-markovian 
approaches, including the numerical renormalisationgroup 
(NRG) U, the numerical hierarchy technique (NHT) JfJ], and 
the numerical path integral (NPI) |7]. There are, however, lim- 
itations concerning the allowed spectral density of the bath 
and these techniques are expected to become less efficient 
with decreasing temperatures and highly structured environ- 



ments. 

Given the general lack of information about the real protein 
spectral densities in PPCs, a technique is required that can 
simulate EET for arbitrary spectral densities and coupling 
strengths, thus allowing experiments carried out under dif- 
ferent conditions, including low temperatures, to be analyzed 
within one framework. A similar situation is encountered in 
solid state qubit implementations, where the microscopic de- 
tails behind the dephasing caused by spurious two level fluc- 
tuators still needs to be fully clarified fl 1 Oil . 
Here we address these issues by developing a numerically ex- 
act and efficient method that combines time-adaptive density 
matrix renormalisation (t-DMRG) with analytical transforma- 
tion techniques based on the theory of orthogonal polynomi- 
als. This method delivers both controllable accuracy and nu- 
merical efficiency, without any restrictions on the complexity 
or strength of the system-environment coupling. Even though 
we will present numerical examples of richly structured en- 
vironments used in the PPC literature, it should be empha- 
sized that this new simulation tool is completely general, and 
can be applied to any system linearly coupled to bosonic or 
fermionic environments. Our method also provides complete 
information about the evolving state of the environment, and 
opens the door to detailed studies of the system-bath correla- 
tions which give rise to long-lasting coherences, entanglement 
and other novel effects. 



The Model- We demonstrate our approach to open-system 
dynamics by considering an elementary model of a PPC, a 
dimer molecule consisting of two pigments (henceforth re- 
ferred to as 'sites'). The dimer's internal dynamics are de- 
scribed by a Hamiltonian Hs = ^kcr-y z + ^o-2 Z +J{cri+o'2- + 
C2+ cr i-)> where , <r,-_ , <7, 2 are standard Pauli creation, an- 
nihilation, and z matrices for the ith site of the dimer. The 
spin-down state a- lz \ \j) = — | li) represents the ground state 
of the site and the spin-up state represents a single local excita- 
tion which can hop between sites with a tunneling amplitude 
J. These sites interact with local independent environments 
modeled as continuous baths of harmonic oscillators, and the 
site-environment interaction Hj and the environment Hamit- 
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FIG. 1 : (a) The standard Hamiltonian has each site of the dimer in- 
teracting with its surrounding bath in a star-like configuration, (b) 
After a unitary transformation of the bosonic modes only, an infi- 
nite harmonic 1 — D chain Hamiltonian is generated with frequen- 
cies e„ and nearest-neighbour couplings t n . The dynamics of the 
dimer is the same in both (a) and (b), but the Hamiltonian of (b) 
can be efficiently simulated using t-DMRG. (c) Intersite population 
transfer rate V as a function of reorganisation energy A for the spec- 
tral density of Eq. (5) at T = OK. Simulation parameters are 
J = 20cm~ 1 ,ei — (.2 — 100cm -1 and 7 = 53cm _1 . Due to 
the presence of multiple timescales in the population transfer, T was 
estimated as the average decay rate of the populations at time 1.3 ps. 



Ionian Hb can be written as 



** = E 



1 + cr,: 



i=l,2 



h{k){ ai {k) + a\(k))dk, (1) 



H B = I 9(k)a\(k)a t (k)dk, 



(2) 



where the environmental modes for each site i are de- 
scribed by creation and annihilation operators a] (k), dj (k) re- 
spectively, and satisfy the continuum commutation relation 
[aj(fc), aj(fc')] = 5ijS(k — k'). For simplicity we assume 
that both baths have identical dispersions g(k) and coupling 
strengths h(k), though our method could easily accommo- 
date different coupling structures on each site. We further 
assume that the spectrum of the bath frequencies is limited 
by a high-frequency cut-off ui c . The form of the coupling 
is chosen so that it is zero when a site is in the ground 
state. The effects of the site-bath interactions are completely 
determined by specifying the spectral function J(u>) of the 



bath lfl7ll . Introducing the inverse function for the disper- 
sion g~ 1 (k), J(ui) is defined for the continuous oscillators by 

J(w) = Trh 2 (g- 1 (u)) d9 ~^}" ) |fTHl . A given J(w) does not 
uniquely specify both g(k) and h{k), but we now show how 
the choice g(k) = oj c k (corresponding to a linear sampling 
of the spectral density) allows the original Hamiltonian to be 
mapped into a ID chain form suitable for efficient DMRG 
simulation. Note that the specific form of the system oper- 
ator in Eq. (Q]) is of no importance for the following trans- 
formations, and the mapping will also apply, for example, for 
environments that induce or suppress spontaneous decay of 
excitations such as those found in photonic crystals IU2I1 . 

We can represent the continuum operators by a new set of 
bosonic operators labelled by a discrete index. These new 
bosonic operators are given by bi n = J Q U n (k)ai(k)dk and 

b\ n = ^ U n {k)a\{k)dk, where n = 0, 1,2.. .00. We now 
choose U n (k) = h(k)p^ l 1 TT n (k), where p n is a normali- 
sation constant, and Tr n (k) is the nth monic polynomial of 
a sequence of orthogonal polynomials defined by the inner 
product Jg 1 h 2 (k)-K n (k)n m (k)dk = p^Snm. With the choice 
g(k) = uj c k, we have h 2 (k) = ■n~ 1 uj c J(g(k)). Using the 
properties of orthogonal polynomials it is easy to show that 
the transformation must be orthogonal and that the new modes 
obey bosonic commutation relation [6„,&JJ = S nm 11411 . Ex- 
pressing H in terms of the new modes, a ID chain Hamilto- 
nian of the form shown in Fig[T]is generated, and is given by 
H = Hs + Hie + He, 

i=l,2 



H. 



CI 



[ (ho + b\ ), 



(3) 



Hr 



1 = 1,2 71 = 



{U n b\ n+1 b in + h.c). (4) 



The new modes of the chain have frequencies ej n and 
nearest-neighbour couplings tj n . The coupling of site j 
to the first member of the chain is determined by 77 = 
J c J{oS) duj. The nearest-neighbour chain structure is a di- 
rect consequence of using the basis of orthogonal polyno- 
mials. When the new operators are substituted into Eq.(Q]) 
the effective coupling strength of the nth mode of the chain 
to the site is oc J Q h 2 {k)-K n {k) dk, which is only non- 
zero for n — due to the orthogonality of the w n (x) 
and the fact that ttq(x) = 1, Similarly, when the new 
modes are substituted into Eq.©, modes n and m couple to 
each other with a strength oc h 2 (k) kn n (k)TT m (k) dk = 

h 2 (k) (a„7r„(fc)+/3„7r«-i(fc)+7r n+ i(fc))7r m (A:) dk, where 
we have used the standard three-term recurrence relation 
obeyed by all monic orthogonal polynomials to replace the 
product kir n (k). After this step, orthogonality guaranties 
nearest-neighbour couplings only, and the whole system of 
sites and oscillators can be arranged in an infinite ID chain 
as shown in Figfl] Note that we have chosen a dimer setting 
only for simplicity, and that systems with more components 
can also be treated efficiently using known generalizations of 
DMRG. 
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This chain mapping is exact and does not require any dis- 
cretisation of the spectral density before being applied.The 
transformation, which also works for fermions, and its analyt- 
ical properties will presented rigorously in Ref. 01411 . where 
analytical formulae for e„, t n can be found for both the linear 
DMRG and the logarithimically discretised NRG chain map- 
pings for J(ui) oc lu s . 1 1 ill . When analytical results are not 
available, the use of orthogonal polynomials also offers sig- 
nificant numerical advantages, as it can be shown that Ci n and 
t i n are simply related to the recurrence coefficients a„ , (3 n in- 
troduced above, and can be determined accurately using very 
stable numerical algorithms such as the ORTHOPOL package 
lfl3l [l4ll . Interestingly, mathematical theorems on these trans- 
formations demonstrate that for any smooth spectral density 
which is strictly bigger than zero in its domain with cut-off 
w c , ti n — > u c /2,ti n — > w c /4 as n — > 00. This result is the 
mathematical expression of the fact that for a smoth spectral 
density we expect that excitations are eventually lost to the 
environment irreversibly. In the linear chain picture this irre- 
versibility arises because a translationally invariant harmonic 
chain with this ratio of e„ to t n has a gapless spectrum and can 
carry away arbitrary excitations with energies up to u> c . Thus 
there is no unlimited build-up of excitations in the chain and 
correlations are expected to be bounded. This in turn suggests 
that DMRG will converge quickly 111611 . 

The overdamped brownian oscillator spectral density - We 
now consider some specific spectral densities of relevance for 
PPCs in photosynthetic organisms. To start with, we look at 
the overdamped Brownian oscillator spectral density which 
has previously been studied by Ishizaki and Fleming in the 
high temperature limit The bath in our simulations is at 
T = K for simplicity, although we can extend our method 
to finite temperatures using standard extensions of DMRG to 
mixed states 111 511 . In our notation, the overdamped brownian 
oscillator spectral density has a simple Ohmic form, 



J(w) 



8A7W 



(5) 



where A is the reorganisation energy of the bath, defined by 
A = J Q " J(u)ui~ 1 dui, and is taken as our measure of the 
site-environment coupling strength. The parameter 7 approx- 
imately sets the dynamical response time of the bath, and the 
following simulations use values of 7 which are smaller than 
the dimer energy scales in order to observe non-markovian ef- 
fects @]. In a PPC an initial excitation is injected suddenly 
onto a site from either another PPC or the capture of a photon. 
The initial condition of the subsequent dynamics at T = K 
is the separable state of zero excitations in the bath, and one 
excitation on site one. This initial condition is used through- 
out this work, although our approach allows us to simulate the 
evolution of any initial state, including highly correlated states 
of the system and bath. 

Figure (Q~|)c shows the average transfer rate of the popula- 
tion initially on site 1 to site 2 as a function of A. We see here 
that the transfer rate is a non-monotonic function of A, with 
an optimum transfer rate occurring for A rj 80cm -1 . Ishizaki 
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FIG. 2: Evolutions of the population on site 1 for the spectral den- 
sity of Eq. |T5} at T = K., and various reorganisation energies A. 
Simulation parameters are J = 100cm" 1 , ei - e 2 = 100cm" 1 and 
53cm" 1 
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and Fleming studied the same dimer and bath parameters at 
300 K using their NHT and found completely incoherent and 
exponential relaxation for A > 5cm" 1 , with the resulting in- 
tersite rate also being a non-monotonic function of A with a 
maximum at A w 10 - 20cm" 1 . This non-monotonicity is a 
signature of the non-perturbative nature of both simulations, 
and in the cases of the NHT, it shows that the method interpo- 
lates correctly between a weak-coupling Redfield-like theory 
and the incoherent Forster theory where multi-phonon effects 
dominate @]. However, the DMRG population evolutions 
(not shown) show that coherence is more robust at T = 
K, with weakly-damped oscillations for A < 30cm" 1 , fol- 
lowed by incoherent relaxation with multiple decay rates for 
30 > A < 50cm" 1 . For larger A the localisation effect dis- 
cussed below comes increasingly into play. 

Figure (0 shows the population on site 1 as a function of 
time for various As. For A < 100cm" 1 we find damped oscil- 
lations which persist for at least 1 ps. For larger A, coherent 
dynamics are always seen for a few hundred femtoseconds 
before the dynamics becomes incoherent, although as A in- 
creases the duration of coherent motion becomes shorter. For 
A > 200cm" 1 the incoherent relaxation rate decreases dra- 
matically, and an increasingly large population is trapped on 
site 1 over the timescale of the simulations. This quantum- 
zeno-like phenomenon may be related to the well-studied lo- 
calisation transition found in Ohmic and sub-Ohmic spin- 
boson models at T = OK 111711 . This is another non- 
perturbative feature of the dynamics, and similar dynamics 
have also recently been observed in NRG and NPI studies of 
the sub-Ohmic spin-boson model J^,[T§]. 

Other spectral densities - We now demonstrate the versa- 
tility of our method w.r.t. the microscopic system-bath inter- 
actions by considering a much more complex and structured 
environmental spectral function taken from a recent study of 
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photosynthetic EET. In Ref lfl9ll . Adolphs and Renger use a 
combination of super-ohmic densities and a coupling to a sin- 
gle effective high-energy mode to model the environment. In 
our notation this spectral function can be written as, 



JM = 



2ttA flOOOw 5 



4.3w 5 



9!(1000o;f 



4.3c 2 5 ) 



+ 4:TtSh^hS(uj — uj h ), 



(6) 



where we have kept the relative contributions of the two 
continuous parts of the spectral density as they are in lfl9ll . 
but have also introduced an overall reorganisation energy A 
to be used as a free parameter. The coupling to the high- 
energy mode is fixed, and the parameters of the simulation 
are J = 100cm -1 , ei — e? = 100cm -1 , wi = 0.5cm -1 



1000cm" 1 and S H = 0.22 



, ei - e 2 

1.95cm- 1 , u H = 180cm -: 

lfl9ll . With these values the continuous part of J(uj) extends 
over a frequency range of about 900cm -1 , and ojh is almost 
resonant with the energy difference ( 224cm -1 ) of the dimer 
eigenstates of H s as the coupling strength of this mode to a 
site is 84cm -1 . 

The chain transformation and DMRG method offers numer- 
ical advantages over some other techniques for spectral func- 
tions which contain delta functions or damped resonances, as 
strong coupling to such modes of the environment do not have 
to be considered part of the quantum system itself. A simple 
modification of the chain parameters allows simulation of an 
arbitrary number of such mode interactions with no increase 
in the complexity of the simulation. Coupling to undamped 
modes with frequencies comparable to or smaller than the 
dimer energies have to be consider as part of the system in 
other approaches like NPI, and this always increases the sim- 
ulation complexity. 

The interaction with the near-resonant oscillator has a pro- 
nounced effect on the population dynamics, and Fig. [3] 
shows how this coupling leads to a coherent beating effect 
which periodically suppresses population oscillations for A < 
300cm -1 . In situations where site 2 might transfer population 
to another system, such a coherent suppression of oscillations 
could lead to an enhancement of EET from the dimer to that 
system. As A increases, the continuous part of the spectral 
density dominates the dynamics and we observe qualitatively 
similar behaviour to the dynamics obtained in Fig. [2] We 
note that the trapping-like dynamics for large A is less severe 
for this super-Ohmic J{uS), although the dynamics are still 
strongly non-Markovian for strong coupling. 

A particularly striking feature of Fig. [3]is that in the regime 
of optimal EET (A ~ 100cm -1 ), the high-energy mode leads 
to low amplitude oscillations which persist for at least 1.5 ps. 
When the high-energy mode is decoupled, coherent oscilla- 
tions vanish for A = 100cm -1 after just 0.3 ps. Experimen- 
tal observation of such persistent undamped oscillations after 
a fast population transfer could thus indicate the presence of 
discrete high-energy modes in the environment of PPCs, and 
could be a useful signature for determining realistic J(u>)s in 
PPCs. 
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FIG. 3: Evolutions of the population on site 1 for the spectral func- 
tion of Eq. l[6} at various reorganization energies A and T = K. 
Dimer parameter are J = 100cm -1 , ei -e 2 = 100cm -1 . Dashed 
line shows how the dynamics when the high-energy mode is decou- 
pled. 



Discussion & Conclusions - Combining an analytical chain 
transformation with numerically exact t-DMRG methods, we 
have introduced an efficient method for the simulation of 
archetypal models of open quantum systems. We have pre- 
sented EET simulations of a dimer molecule in the presence 
of a variety of environmental coupling structures and observed 
a rich array of non-perturbative and non-markovian effects in 
the dynamics. Having demonstrated the viability of the tech- 
nique, we plan to explore the important question of the role 
and persistence of coherence in biological EET, using existing 
DMRG techniques to extend our method to both finite temper- 
atures, multiple sites, and spatially-correlated baths. As men- 
tioned previously, the method is general and we expect it to 
be relevant for a variety of problems in condensed matter and 
quantum information, including fundamental studies of deco- 
herence, such as the spin-boson model, and the boundary of 
quantum and classical physics. 
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